Chapter 4 



Local stability 



(R. Mainieri and P. Cvitanovic) 

So far we have concentrated on describing the trajectory of a single initial 
point. Our next task is to define and determine the size of a neighborhood 
of x{t). We shall do this by assuming that the flow is locally smooth and by 
describing the local geometry of the neighborhood by studying the flow linearized 
around x(t). Nearby points aligned along the stable (contracting) directions remain 
in the neighborhood of the trajectory x(t) = f'(xo); the ones to keep an eye on are 
the points which leave the neighborhood along the unstable directions. As we shall 
demonstrate in chapter 18, the expanding directions matter in hyperbolic systems. 
The repercussions are far-reaching. As long as the number of unstable directions 
is finite, the same theory applies to finite-dimensional ODEs, state space volume 
preserving Hamiltonian flows, and dissipative, volume contracting infinite-dim- 
ensional PDEs. 

In order to streamline the exposition, in this chapter all examples are collected 
in sect. 4.8; you can get to them and back to the text by clicking on the [example] 
links, such as 

example 4.8 
p. 90 



4.1 Flows transport neighborhoods 

As a swarm of representative points moves along, it carries along and distorts 
neighborhoods. The deformation of an infinitesimal neighborhood is best un- 
derstood by considering a trajectory originating near xq = x(0), with an initial 
infinitesimal displacement 6x(0). The flow then transports the displacement 8x{t) 
along the trajectory x(xq, t) - f'(xo). 
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4.1.1 Instantaneous shear 

The system of linear equations of variations for the displacement of the infmites- 
imally close neighbor x + 6x follows from the flow equations (2.7) by Taylor 
expanding to linear order 

Zdvj 
- — 8x; . 
OXj 

J 1 

The infinitesimal displacement 8x is thus transported along the trajectory x(x , t), 
with time variation given by 



-5x i (x ,t) = ± j —(x) 

j 1 



6xj(x ,t). (4.1) 

x=x(x ,t) 



As both the displacement and the trajectory depend on the initial point xq and the 
time t, we shall often abbreviate the notation to x(xo, t) —> x(t) — > x, 6xj(xo, t) — > 
Sxi(t) —> 8x in what follows. Taken together, the set of equations 

Xi = Vj(x) , 5xj = 2^ Ajj(x)Sxj (4.2) 
j 



governs the dynamics in the tangent bundle (x, 5x) e TM obtained by adjoining 
the (/-dimensional tangent space 5x e TM X to every point x € At in the d-dim- 
ensional state space M c W 1 . The stability matrix (velocity gradients matrix) 



A u (x) = ^ (4.3) 

OXj 

describes the instantaneous rate of shearing of the infinitesimal neighborhood of 
x(t) by the flow. 



example 4.1 
p. 90 



4.1.2 Finite time linearized flow 

By Taylor expanding a. finite time flow to linear order, 



fjixo + Sx) = f[(x ) + V — dxj + ■■■ , (4.4) 

OXq : 

J J 
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Figure 4.1: A local frame is transported along the or- r §b]0 
bit and deformed by Jacobian matrix J'. As J' is not 
self-adjoint, initial orthogonal frame is mapped into a 
non-orthogonal one. 

rgb]0,0,0x(0) 

one finds that the linearized neighborhood is transported by the Jacobian matrix remark 4.1 




t t ® x i 

6x(t) = J (x )dx , J u (x ) - - — 

1 dx Qj 



, J°(x Q ) = 1 . 



(4.5) 



x=x(t) 



For example, in 2 dimensions the Jacobian matrix for change from initial to final 
coordinates is 



ar \ ( dx \ 

J - (\X) - — Q y Sy 

d(x ,yo) {or Wo J 



The map f is assumed invertible and differentiable so that J' exists. For suf- 
ficiently short times J' remains close to 1, so det/ f > 0. By continuity det/ r 
remains positive for all times t. However, for discrete time maps, det J' can have 
either sign. 



4.1.3 Covariant frames 



/ describes the deformation of an infinitesimal neighborhood at finite time t in the 
co-moving frame of x(t), or, in the language of sect. 2.2.1, the Jacobian matrix 
maps the initial Lagrangian coordinate frame into the current Eulerian coordinate 
frame. This deformation of an initial frame at xq into a non-orthogonal frame at 
x(t) is described by the eigenvectors and eigenvalues of the Jacobian matrix of the 
linearized flow (see figure 4.1), 

fe U) = Aje {j) , j=l,2,---,d. (4.6) 

Throughout this text the symbol will always denote the &th eigenvalue (some- 
times referred to as the multiplier) of the finite time Jacobian matrix J'. Symbol 
A (k) will be reserved for the &th stability exponent, with real part //^ and phase 

A, = e aik) A {k) = v& + iJ® . (4.7) 



As J' is a real matrix, its eigenvalues are either real or come in complex pairs, 



{A,,A, +1 }-{ e '^ ft,+ - W ),e"^ , -- (A ' ) >} 
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with magnitude |A^| and phase a>. The phase describes the rotation velocity 
in the plane defined by the corresponding pair of real orthogonal eigenvectors, 
{Ree (/:) , Ime^}, with one period of rotation given by T = 2nj<x> . 



example 4.4 
p. 91 

J'{xq) depends on the initial point xq and the elapsed time t. For notational 
brevity we omitted this dependence, but in general both the eigenvalues and the 
eigenvectors, Aj = Aj(xq, t) , • • • , e' 7 - 1 = e^(jco, t) , also depend on the trajectory 
traversed and the choice of coordinates. We shall refer to the eigenvectors as 
covariant vectors. 

Nearby trajectories separate along the unstable directions, approach each other 
along the stable directions, and change their distance along the marginal direc- 
tions at a rate slower than exponential, corresponding to the eigenvalues of the Ja- 
cobian matrix with magnitude larger than, smaller than, or equal to 1 . In the litera- 
ture, the adjectives neutral, indifferent, center are often used instead of 'marginal'; 
(attracting) stable directions are sometimes called 'asymptotically stable,' and so 
on. 

One of the preferred directions is what one might expect, the direction of the 
flow itself. To see that, consider two initial points along a trajectory separated 
by infinitesimal flight time St: 6xq = f 6t (xo) - xq = v(xo)5t . By the semigroup 
property of the flow, f t+6t = f 6t+t , where 

XSt+t 
drv(x(T)) + f\xQ) = Stv(x(t)) + f(x ) . 

Expanding both sides of f'{f 5t {xQ)) = f St (f'(xo)), keeping the leading term in 
St, and using the definition of the Jacobian matrix (4.5), we observe that J'(xq) 
transports the velocity vector at xq to the velocity vector at x{t) (see figure 4.1): 

v(x(t)) = J ! (x )v(x ). (4.8) 




4.2 Computing the Jacobian matrix 

As we started by assuming that we know the equations of motion, from (4.3) we 
also know stability matrix A, the instantaneous rate of shear of an infinitesimal 
neighborhood 8xj{t) of the trajectory x{t). What we do not know is the finite time 
deformation (4.5), so our next task is to relate the stability matrix A to Jacobian 
matrix J* . On the level of differential equations the relation follows by taking the 
time derivative of (4.5) and replacing 6x by (4.2) 

d dJ* , 

— 8x(t) - —— 6xq - A Sx(t) = AJ Sx . 

dt dt 



stability - 17apr2013 



ChaosBook.org versionl4.4.1, Apr 21 2013 



CHAPTER 4. LOCAL STABILITY 



81 



Hence the d 2 matrix elements of Jacobian matrix satisfy 'tangent linear equations,' 
the linearized equations (4.1) 

—f(x ) = A(x) J'(x ) , x = f(x ) , initial condition J°(x ) = 1 ■ (4-9) 
at 

For autonomous flows, the matrix of velocity gradients A(x) depends only on x, 
not time, while J' depends on both the state space position and time. Given a 
numerical routine for integrating the equations of motion, evaluation of the Jaco- 
bian matrix requires minimal additional programming effort; one simply extends 
the cf-dimensional integration routine and integrates concurrently with f'(xo) the 
d 2 elements of J'{x$). The qualifier 'simply' is perhaps too glib. Integration will 
work for short finite times, but for exponentially unstable flows one quickly runs 
into numerical over- and/or underflow problems. For high-dimensional flows the 
analytical expressions for elements of A might be so large that it fits on no com- 
puter. Further thought will have to go into implementation this calculation. chapter 26 

So now we know how to compute Jacobian matrix /' given the stability matrix 
A, at least when the d 2 extra equations are not too expensive to compute. Mission 
accomplished. 




fast track: 
chapter 7, p. 125 



And yet... there are mopping up operations left to do. We persist until we de- 
rive the integral formula (4.19) for the Jacobian matrix, an analogue of the finite- 
time 'Green's function' or 'path integral' solutions of other linear problems. 

We are interested in smooth, differentiable flows. If a flow is smooth, in a suf- 
ficiently small neighborhood it is essentially linear. Hence the next section, which 
might seem an embarrassment (what is a section on linear flows doing in a book 
on nonlinear dynamics?), offers a firm stepping stone on the way to understanding 
nonlinear flows. Linear charts are the key tool of differential geometry, general 
relativity, etc., so we are in good company. If you know your eigenvalues and 
eigenvectors, you may prefer to fast forward here. 

fast track: 
sect. 4.4, p. 84 

4.3 A linear diversion 

Linear is good, nonlinear is bad. 
— Jean Bellissard 

Linear fields are the simplest vector fields, described by linear differential equa- 
tions which can be solved explicitly, with solutions that are good for all times. 
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The state space for linear differential equations is M = W 1 , and the equations of 
motion (2.7) are written in terms of a vector x and a constant stability matrix A as 

x = v{x)=Ax. (4.10) 

Solving this equation means finding the state space trajectory 

x{t) = (x l (t),x 2 (t),...,x d (t)) 

passing through a given initial point xq. If x(t) is a solution with x(0) = xq and 
y(t) another solution with y(0) = y , then the linear combination ax(t) + by(t) with 
a, b € R is also a solution, but now starting at the point ax^ + by^. At any instant 
in time, the space of solutions is a <i-dimensional vector space, spanned by a basis 
of d linearly independent solutions. 

How do we solve the linear differential equation (4.10)? If instead of a matrix 
equation we have a scalar one, x - Ax, the solution is x(t) = e tA xo . In order 
to solve the cf-dimensional matrix case, it is helpful to rederive this solution by 
studying what happens for a short time step St. If at time t = the position is x(0), 
then 

x{5t) - x(0) 

= Ax(0) , (4.11) 

St 

which we iterate m times to obtain Euler's formula for compounding interest 

(t V" 
1 + -A) x(0) « e a x(0) . (4.12) 
m I 

The term in parentheses acts on the initial condition x(0) and evolves it to x(t) by 
taking m small time steps 8t = t/m. As m —> oo, the term in parentheses converges 
to e a . Consider now the matrix version of equation (4.11): 



x(6t) - x(0) 

KJ =Ax(0). (4.13) 

ot 



A representative point x is now a vector in R acted on by the matrix A, as in 
(4. 10). Denoting by 1 the identity matrix, and repeating the steps (4.1 1) and (4. 12) 
we obtain Euler's formula for the exponential of a matrix: 

(t \ m 
1 + -A . (4.14) 
m I 

We will find this definition the exponential of a matrix helpful in the general case, 
where the matrix A = A(x(t)) varies along a trajectory. 
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As the system is a linear, a superposition of any two solutions to x(t) = J'x(0) 
is also a solution. One can take any d independent initial states, x (1) (0), x®(0), 
. . . , x (d \0), assemble them as columns of a matrix (D(0), and formally write the 
solution for an arbitrary initial condition projected onto this basis, 

x(t) = Oa)O(0) _1 x(0) (4.15) 

where <£(?) - [x (1 \t), x (2) (f), • • • , x (d) (t)\. Q>(t) is called the fundamental matrix of 
the system, and the Jacobian matrix /' = OtOOCO) -1 can thus be fashioned out of 
d trajectories {x^\t)}. Numerically this works for sufficiently short times. 

Now that we have some feeling for the qualitative behavior of eigenvectors and 
eigenvalues of linear flows, we are ready to return to the nonlinear case. How do 
we compute the exponential (4. 14)? 

example 4.2 TKly^* fast track: section 5.2.2 

p. 90 fc^P S ect. 4.4, p. 84 

Henriette Roux: So, computing eigenvalues and eigenvectors seems like a good 
thing. But how do you really do it? 

A: Any text on numerics of matrices discusses how this is done; the keywords are 
Gram-Schmidt, and for high-dimensional flows 'Krylov subspace' and 'Arnoldi 
iteration' . Conceptually (but not for numerical purposes) we like the economical 
description of neighborhoods of equilibria and periodic orbits afforded by projec- 
tion operators. The requisite linear algebra is standard, but usually not phrased in 
language of projection operators. As this is a bit of sidetrack that you will find 
confusing at the first go, it is relegated to appendix C. 

Henriette Roux: In my 61,506-dimensional computation of a Navier-Stokes equi- 
librium I generated about 30 eigenvectors before I wanted to move on. How many 
of these eigenvectors are worth generating for a particular solution and why? chapter ?? 

A: A rule of the thumb is that you need all equilibrium eigenvalues / periodic orbit 
Floquet exponents with positive real parts, and at least those negative exponents 
whose magnitude is less or comparable to the largest expanding eigenvalue. More 
precisely; keep adding the next least contracting eigenvalue to the sum of the 
preceding ones as long as the sum is positive (Kaplan- Yorke criterion). Then, just 
to be conservative, double the number of eigenvalues you keep. You do not need to 
worry about the remaining (60 thousand!) eigen-directions for which the negative 
eigenvalues are of larger magnitude, as they always contract: nonlinear terms 
cannot mix them up in such a way that expansion in some directions overwhelms 
the strongly contracting ones. 
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4.4 Stability of flows 

How do you determine the eigenvalues of the finite time local deformation J' for 
a general nonlinear smooth flow? The Jacobian matrix is computed by integrating 
the equations of variations (4.2) 



x(t) - f{xo), 6x(xo,t) = J\xq) 6x(xq, 0), 



(4.16) 



The equations are linear, so we should be able to integrate them-but in order to 
make sense of the answer, we derive this integral step by step. 

Consider the case of a general, non-stationary trajectory x(t). The exponential 
of a constant matrix can be defined either by its Taylor series expansion, or in 
terms of the Euler limit (4.14): 



JA 



yL A k = lim (i + La) . 

k=0 



(4.17) 



Taylor expanding is fine if A is a constant matrix. However, only the second, 
tax-accountant's discrete step definition of an exponential is appropriate for the 
task at hand, as for a dynamical system the local rate of neighborhood distortion 
A{x) depends on where we are along the trajectory. The linearized neighborhood 
is multiplicatively deformed along the flow, and the m discrete time-step approx- 
imation to f is therefore given by a generalization of the Euler product (4.17): 



/ = lim li (1 + 5tA{x n )) = lim FT e StA(x " 

m—tco 1 1 m—>oo 1 1 



(4.18) 



= lim e' 

m—yco 



StA(x m ) e StA(x m -i) 



, e StA(x 2 ) e 6tA(x 1 ) 



where St — (t - to)/m, and x n - x(to + nSt). Indexing of the products indicates that 
the successive infinitesimal deformation are applied by multiplying from the left. 
The m — > oo limit of this procedure is the formal integral 



4-(*o) = 



(4.19) 



where T stands for time-ordered integration, defined as the continuum limit of 
the successive left multiplications (4.18). This integral formula for J' is the exercise 4.5 
main conceptual result of this chapter, the finite time companion of the differential 
definition (4.9). The definition makes evident important properties of Jacobian 
matrices, such as that they are multiplicative along the flow, 



J t+t, (x) = j\x') J\x), where x' = f(x ) , 



(4.20) 
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an immediate consequence of time-ordered product structure of (4.18). However, 
in practice / is evaluated by integrating (4.9) along with the ODEs that define a 
particular flow. 



4.5 Stability of maps 



The transformation of an infinitesimal neighborhood of a trajectory under the iter- 
ation of a map follows from Taylor expanding the iterated mapping at finite time 
n to linear order, as in (4.4). The linearized neighborhood is transported by the 
Jacobian matrix evaluated at a discrete set of times n - 1,2, . . ., 



(4.21) 



x=x 



As in the finite time case (4.7), we denote by A* the kth eigenvalue or multiplier 
of the finite time Jacobian matrix J" . There is really no difference from the con- 
tinuous time case, other than that now the Jacobian matrix is evaluated at integer 
times. 



example 4.9 
p. 95 

The formula for the linearization of nth iterate of a J-dimensional map 

J'\x ) = /(*„_!) • • • J(xOJ(xo) , xj = f j (x Q ) , (4.22) 

in terms of single time steps Jj; = dfj/dxi follows from the chain rule for func- 
tional composition, 



d ^ dfj(y) 

dxi frf dy k 



dfk(x) 

y=f(x) dx < 



If you prefer to think of a discrete time dynamics as a sequence of Poincare sec- 
tion returns, then (4.22) follows from (4.20): Jacobian matrices are multiplicative 
along the flow. 



example 4.10 7^-^-= fast track: 

p. 96 chapter 7, p. 125 
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Figure 4.2: If x(t) intersects the Poincare section 
P at time t, the nearby x(t) + 6x(t) trajectory inter- 
sects it time r + St later. As ((/' ■ v'dt) = -((/'• 
J dx), the difference in arrival times is given by St = 
-(U' ■JSx)/(U' -v'). 




4.6 Stability of Poincare return maps 



(R. Paskauskas and P. Cvitanovic) 



We now relate the linear stability of the Poincare return map P : P — > P denned 
in sect. 3.1 to the stability of the continuous time flow in the full state space. 

The hypersurface P can be specified implicitly through a function U (x) that is 
zero whenever a point x is on the Poincare section. A nearby point x + 8x is in the 
hypersurface P if U(x+5x) = 0, and the same is true for variations around the first 
return point x' = x(t), so expanding U{x') to linear order in variation 8x restricted 
to the Poincare section, and applying the chain rule leads to the condition 



In what follows Uj - djU is the gradient of U defined in (3.3), unprimed quantities 
refer to the starting point x = xq € P, v = v(xrj), and the primed quantities to the 
first return: x' = x(t), V = v(x'), U' = U(x'). For brevity we shall also denote 
the full state space Jacobian matrix at the first return by J = J t (xq). Both the first 
return x' and the time of flight to the next Poincare section t(x) depend on the 
starting point x, so the Jacobian matrix 



with both initial and the final variation constrained to the Poincare section hyper- 
surface P is related to the continuous flow Jacobian matrix by 



The return time variation dr/dx, figure 4.2, is eliminated by substituting this ex- 
pression into the constraint (4.23), 




(4.23) 




(4.24) 




= diU' Ju + (/ • dU')— , 
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yielding the projection of the full space J-dimensional Jacobian matrix to the 
Poincare map (d- l)-dimensional Jacobian matrix: 



Jii = \8ik 



V' : dlrW 



7\\ Jk j 



(4.25) 



(v' -dU'), 

Substituting (4.8) we verify that the initial velocity v(x) is a zero-eigenvector of J 
Jv = , (4.26) 



so the Poincare section eliminates variations parallel to v, and J is a rank (d- 1)- 
dimensional matrix, i.e., one less than the dimension of the continuous time flow. 



4.7 Neighborhood volume 

Consider a small state space volume AV = d x centered around the point xo at 
time t = 0. The volume A V around the point xf = x{t) time t later is 



secuori d.^ 
remark 6.1 



, AV 

AY = AV 

AV 



, dx' 




det- 


AV - 


ox 



(4.27) 



so the |det /| is the ratio of the initial and the final volumes. The determinant 

det/'(xo) = nf=i A,(xo, t) is the product of the Jacobian matrix eigenvalues. We 

shall refer to this determinant as the Jacobian of the flow. The Jacobian is easily exercise 4.1 

evaluated. Take the time derivative, use the / evolution equation (4.9) and the 

matrix identity In det J = tr In J: 

d d d 1 . 

— In AV(t) = — In det J = tr — In J = tr -J = tr A = dm . 

dt dt dt J 



(Here, as elsewhere in this book, a repeated index implies summation.) Integrate 
both sides to obtain the time evolution of an infinitesimal volume 



det J'(xq) = exp 



f 

Jo 



<frtrA(x(T)) 



f 

Jo 



= exp dTdiVi(x(T)) 



(4.28) 



As the divergence 3;V, is a scalar quantity, the integral in the exponent (4. 19) needs 
no time ordering. So all we need to do is evaluate the time average 



Y rt d 

dm = lim - I dr V A u {x{t)) 
t Jo 



t 



;=1 



~\Ai{x Q ,t) =Y J A (i \xo,t) 



(4.29) 



;=1 
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along the trajectory. If the flow is not singular (for example, the trajectory does 
not run head-on into the Coulomb 1/r singularity), the stability matrix elements 
are bounded everywhere, |Aj y -| < M , and so is the trace 2,-A/,\ The time integral 
in (4.29) thus grows at most linearly with t, d(Vi is bounded for all times, and 
numerical estimates of the t — > oo limit in (4.29) are not marred by any blowups. 



Even if we were to insist on extracting 3,V; from (4.18) by first multiplying 
Jacobian matrices along the flow, and then taking the logarithm, we can avoid ex- 
ponential blowups in J f by using the multiplicative structure (4.20), det J (jco) = 
det/' (jc') det J\xq) to restart with J°(x') = 1 whenever the eigenvalues of J\xq) 
start getting out of hand. In numerical evaluations of stability exponents the sum section 6.2 
rule (4.29) can serve as a helpful check on the accuracy of the computation. 

The divergence 3,-V; characterizes the behavior of a state space volume in the 
infinitesimal neighborhood of the trajectory. If 3/V; < 0, the flow is locally con- 
tracting, and the trajectory might be falling into an attractor. If 3,-Vf(ac) < , for 
all x € At, the flow is globally contracting, and the dimension of the attractor is 
necessarily smaller than the dimension of state space At. If 3/V; = 0, the flow 
preserves state space volume and det J* — 1. A flow with this property is called 
incompressible. An important class of such flows are the Hamiltonian flows 
considered in sect. 7.3. 

But before we can get to that, Henriette Roux, the perfect student and always 
alert, pipes up. She does not like our definition of the Jacobian matrix in terms of 
the time-ordered exponential (4.19). Depending on the signs of multipliers, the 
left hand side of (4.28) can be either positive or negative. But the right hand side 
is an exponential of a real number, and that can only be positive. What gives? As 
we shall see much later on in this text, in discussion of topological indices arising 
in semiclassical quantization, this is not at all a dumb question. 



Resume 

A neighborhood of a trajectory deforms as it is transported by a flow. In the 
linear approximation, the stability matrix A describes the shearing / compression 
/ expansion of an infinitesimal neighborhood in an infinitesimal time step. The 
deformation after a finite time t is described by the Jacobian matrix 



example 4.8 
p. 95 



where T stands for the time-ordered integration, defined multiplicatively along 
the trajectory. For discrete time maps this is multiplication by time-step Jacobian 
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matrix / along the n points xq, x\, x%, . . ., x n -\ on the trajectory of xq, 

J n (x ) - J(x n -\)J(x n -2) ■ ■ ■ J(xi)J(x ) , 

with J(x) the single discrete time-step Jacobian matrix. In ChaosBook the stability 
multiplier At denotes the kth eigenvalue of the finite time Jacobian matrix J'(xo), 
f/® the real part of kth stability exponent, and co^ its phase, 

A - e t(M+ico ^ 

For complex eigenvalue pairs the 'angular velocity' a> describes rotational motion 
in the plane spanned by the real and imaginary parts of the corresponding pair of 
complex eigenvectors. 

The eigenvalues and eigen-directions of the Jacobian matrix describe the de- 
formation of an initial infinitesimal cloud of neighboring trajectories into a dis- 
torted cloud a finite time t later. Nearby trajectories separate exponentially along 
unstable eigen-directions, approach each other along stable directions, and change 
slowly (algebraically) their distance along marginal, or center directions. The Ja- 
cobian matrix J' is in general neither symmetric, nor diagonalizable by a rotation, 
nor do its (left or right) eigenvectors define an orthonormal coordinate frame. 
Furthermore, although the Jacobian matrices are multiplicative along the flow, in 
dimensions higher than one their eigenvalues in general are not. This lack of 
multiplicativity has important repercussions for both classical and quantum dy- 
namics. 



Commentary 

Remark 4.1 Linear flows. The subject of linear algebra generates innumerable tomes 
of its own; in sect. 4.3 we only sketch, and in appendix C recapitulate a few facts that 
our narrative relies on: a useful reference book is [4.1]. The basic facts are presented at 
length in many textbooks. Frequently cited linear algebra references are Golub and Van 
Loan [4.2], Coleman and Van Loan [4.3], and Watkins [4.4, 4.5]. The standard references 
that exhaustively enumerate and explain all possible cases are Hirsch and Smale [4.6] and 
Arnol'd [4.7]. A quick overview is given by Izhikevich [4.8]; for different notions of orbit 
stability see Holmes and Shea-Brown [4.9]. For ChaosBook purposes, we enjoyed the 
discussion in chapter 2 Meiss [4.10], chapter 1 of Perko [4.11] and chapters 3 and 5 of 
Glendinning [4.12], and liked the discussion of norms, least square problems, and differ- 
ences between singular value and eigenvalue decompositions in Trefethen and Bau [4. 13]. 
Truesdell [2.2] and Gurtin [2.3] are excellent references for the continuum mechanics per- 
spective on state space dynamics. We enjoyed Christov et al. [2.1] easy introduction to 
parallels between dynamical systems and continuum mechanics. 

The nomenclature tends to be a bit confusing. Jacobian matrix (4.5) is sometimes 
referred to as the fundamental solution matrix or simply fundamental matrix, a name in- 
herited from the theory of linear ODEs, or the Frechet derivative of the nonlinear mapping 
f'(x), or the 'tangent linear propagator', or even, by Lorenz [6.2], the 'error matrix'. As 
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/ describes the deformation of an infinitesimal neighborhood at finite time t in the co- 
moving frame of x(t), in continuum mechanics it is called the deformation gradient, or 
transplacement gradient. It is often denoted Df, but for our needs (we shall have to sort 
through a plethora of related Jacobian matrices) matrix notation / is more economical. 
Single discrete time-step Jacobian Jji = dfj/dxi in (4.22) is referred to as the 'tangent 
map' by Skokos [4.16, 4.17]. 

In referring to the 'velocity gradients matrix' or 'velocity gradient tensor' A defined 
in (4.3) as the 'stability matrix' we follow Tabor [4.15]. It is the natural object for study 
of stability of equilibria, time-invariant point in state space; stability of trajectories is 
described by Jacobian matrices. Goldhirsch, Sulem, and Orszag [4.18] call it the 'Hes- 
senberg matrix', and to the equations of variations (4.1) as 'stability equations.' Manos et 
al. [4.19] refer to (4.1) as the 'variational equations'. 

Sometimes A, which describes the instantaneous shear of the neighborhood of x(xq, f), 
is referred to as the 'Jacobian matrix,' a particularly unfortunate usage when one consid- 
ers linearized stability of an equilibrium point (5.1). A is not a Jacobian matrix, just as 
a generator of SO(2) rotation is not a rotation; A is a generator of an infinitesimal time 
step deformation, J St ^ 1 + A5t . What Jacobi had in mind in his 1841 fundamental pa- 
per [4.20] on the determinants today known as 'jacobians' were transformations between 
different coordinate frames. These are dimensionless quantities, while dimensionally Ay 
is l/[time]. 

More unfortunate still is referring to the Jacobian matrix J' = exp(fA) as an 'evolution 
operator,' which here (see sect. 17.2) refers to something altogether different. In this book 
Jacobian matrix /' always refers to (4.5), the linearized deformation after a finite time t, 
either for a continuous time flow, or a discrete time mapping. 



4.8 Examples 

The reader is urged to study the examples collected here. To return back to the 
main text, click on [click to return] pointer on the margin. 



Example 4.1 Rossler and Lorenz flows, linearized: (continued from example 3.5) For 
the Rossler (2.18) and Lorenz (2.13) flows the stability matrices are, respectively 





( 


-1 -1 > 




' —<T cr 
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Arqss — 
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> ^Lor — 


P - z -1 


X 
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x - c , 




, y x 


-b J 



(continued in example 4.5) click to return: p. ?? 



Example 4.2 Jacobian matrix eigenvalues, diagonalizable case: Should we be 
so lucky that A = A D happens to be a diagonal matrix with eigenvalues (A (l \ A (2) , A {d) ), 
the exponential is simply 



f = e tA ° 



( e aW ••• 
■ 



(4.31) 
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Figure 4.3: Streamlines for several typical 2- 
dimensional flows: saddle (hyperbolic), in node (at- 
tracting), center (elliptic), in spiral. 




Next, suppose that A is diagonalizable and that U is a nonsingular matrix that brings it 
to a diagonal formA D = U~ l AU. Then J can also be brought to a diagonal form (insert 
factors 1 = UU~ l between the terms of the product (4. 14)): exercise 4.2 

f = e tA = Ue tAD lT 1 . (4.32) 

The action of both A and J is very simple; the axes of orthogonal coordinate system 
where A is diagonal are also the eigen-directions of J', and under the flow the neigh- 
borhood is deformed by a multiplication by an eigenvalue factor for each coordinate 
axis. 

We recapitulate the basic facts of linear algebra in appendix C. The following 
2-dimensional example serves well to highlight the most important types of linear 
flows: 

Example 4.3 Linear stability of 2-dimensional flows: For a 2-dimensional flow the 
eigenvalues A (l \A {2) of A are either real, leading to a linear motion along their eigen- 
vectors, Xj{t) = xj(Q) exp(f/l ( -' ) ), or a form a complex conjugate pair A (l) = fi + ico, A (2) = 
fx-ico, leading to a circular or spiral motion in the [xi,X2] plane. 

These two possibilities are refined further into sub-cases depending on the 
signs of the real part. In the case of real A {1) > 0, A (2) < 0, x\ grows exponentially 
with time, and x 2 contracts exponentially. This behavior, called a saddle, is sketched in 
figure 4.3, as are the remaining possibilities: in/out nodes, inward/outward spirals, and 
the center. The magnitude of out-spiral \x(t)\ diverges exponentially when /i > 0, and 
in-spiral contracts into (0, 0) when the n < 0, whereas the phase velocity u controls its 
oscillations. 

If eigenvalues A m = A {2) = A are degenerate, the matrix might have two linearly 
independent eigenvectors, or only one eigenvector. We distinguish two cases: (a) 
A can be brought to diagonal form, (b) A can be brought to Jordan form, which (in 
dimension 2 or higher) has zeros everywhere except for the repeating eigenvalues on 
the diagonal, and some 1 's directly above it. For every such Jordan [d xd a ] block there 
is only one eigenvector per block. 

We sketch the full set of possibilities in figures 4.3 and 4.4, and work out in 
detail the most important cases in appendix C, example C.3. click to return 



Example 4.4 In-out spirals. Consider an equilibrium whose stability exponents 
{A W ,A (2) } = \n + i(D,fi - icj] form a complex conjugate pair. The corresponding com- 
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saddle out node in node 



Figure 4.4: Qualitatively distinct types of expo- 
nents , /l <2) ) of a [2x2] Jacobian matrix. 




plex eigenvectors can be replaced by their real and imaginary parts, {e (1) ,e {2) } 
{Re e (1) , Im e (1) }. The 2-dimensional real representation, 



consists of the identity and the generator of SO(2) rotations in the {Ree (1) ,Ime (1) } plane. 
Trajectories x(t) = J' x(0), where (omitting e (3) , e (4) , • • ■ eigen-directions) 

/ = e V = e W cos ^ -**o*\ (433) 
\ sin cot cos cot j 

spiral in/out around (x,y) = (0,0), see figure 4.3, with the rotation period T, and con- 
traction/expansion radially by the multiplier h mdial , and by the multiplier Aj along the 
e 0) eigen-direction per a turn of the spiral: exercise C.1 

T = 2tt/co, A radial = e T f, Aj = e T " m . (4.34) 

We learn that the typical turnover time scale in the neighborhood of the equilibrium 
(x,y) = (0,0) is of order * T (and not, let us say, 1000 T, or 10- 2 rj. A, multipliers 
give us estimates of strange-set thickness in eigen-directions transverse to the rotation 
plane. 



Example 4.5 Stability of equilibria of the Rossler flow. (continued from ex- 
ample 4.1) The Rosier system (2.18) has two equilibrium points (2.19), thesmteise 4.4 
equilibrium (x-,y-,z-), and the outer equilibrium point (x + ,y + ,z + ). Together wittextomise 2.8 
exponents (eigenvalues of the stability matrix), the two equilibria yield quite detailed 
information about the flow. Figure 4.5 shows two trajectories which start in the neigh- 
borhood of the outer '+ ' equilibrium. Trajectories to the right of the equilibrium point '+ ' 
escape, and those to the left spiral toward the inner equilibrium point '- ', where they 
seem to wander chaotically for all times. The stable manifold of outer equilibrium point 
thus serves as the attraction basin boundary. Consider now the numerical values for 
eigenvalues of the two equilibria 

(^/'V- 2 ' + ito {2) ) = (-5.686, 0.0970 + /0.9951 ) 

Gu^./zf + icof ) = (0.1929, -4.596 x lO" 6 + Z5.428) 

Outer equilibrium: The fi+ + i co t2) complex eigenvalue pair implies that neighborhood 
of the outer equilibrium point rotates with angular period T + » |2^/w ( + 2) | = 1.1575. The 
multiplier by which a trajectory that starts near the '+' equilibrium point contracts in the 
stable manifold plane is the excruciatingly slow multiplier A+ « exp( / u <2) r + ) = 0.9999947 
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Figure 4.5: Two trajectories of the Rossler flow initi- 
ated in the neighborhood of the '+' or 'outer' equilib- 
rium point (2.19). (R. Paskauskas) 




per rotation. For each period the point of the stable manifold moves away along the 
unstable eigen-direction by factor A| « expOu^TV) = 1.2497. Hence the slow spiraling 
on both sides of the '+ ' equilibrium point. 

Inner equilibrium: The + iu>_ complex eigenvalue pair tells us that neighbor- 
hood of the '-' equilibrium point rotates with angular period T » l^/to?^ = 6.313, 
slightly faster than the harmonic oscillator estimate in (2.15). The multiplier by which 
a trajectory that starts near the '-' equilibrium point spirals away per one rotation is 
^■radial ~ expfju^T-) = 1.84. The fi_ eigenvalue is essentially the z expansion cor- 
recting parameter c introduced in (2. 1 7). For each Poincare section return, the trajec- 
tory is contracted into the stable manifold by the amazing factor of Ai ~ expQrjT-) = 
lfr 15 - 6 (I). 

Suppose you start with a 1 mm interval pointing in the Ai eigen-direction. After 
one Poincare return the interval is of order of 10~ 4 fermi, the furthest we will get into sub- 
nuclear structure in this book. Of course, from the mathematical point of view, the flow 
is reversible, and the Poincare return map is invertible. (continued in example 1 1.3) 

(R. Paskauskas) 



Example 4.6 Stability of Lorenz flow equilibria: (continued from example 4.1) A 
glance at figure 3.4 suggests that the flow is organized by its 3 equilibria, so lets have 
a closer look at their stable/unstable manifolds. 

The EQo equilibrium stability matrix (4.30) evaluated at x EQo = (0, 0, 0) is block- 
diagonal. The z-axis is an eigenvector with a contracting eigenvalue A (2) = -b. r§naark 9.14 
(4.41) it follows that all [x,y] areas shrink at rate -{a + 1). Indeed, the [x,y] submatrix 

7 -i i < 4 - 36 ) 



has a real expanding/contracting eigenvalue pairA (l ' 3) = -(<r+ 1)/2+ ^(<r- l) 2 /4 + per, 
with the right eigenvectors e (1) , e (3) in the [x,y] plane, given by (either) column of the 
projection operator 

A- - A u) l 1 / -o- - o- \ 

P,- = — -- = -- , ,ffl . »*/e{l,3}. (4.37) 

EQ\,i equilibria have no symmetry, so their eigenvalues are given by the roots 
of a cubic equation, the secular determinant det (A - AI) = 0: 

A 3 + A 2 (cr + b + 1) + Ab(cr + p) + 2crb(p - 1) = . (4.38) 

For p > 24.74, £gi,2 have one stable real eigenvalue and one unstable complex con- 
jugate pair, leading to a spiral-out instability and the strange attractor depicted in fig- 
ure 2.5. 
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Figure 4.6: (a) A perspective view of the lin- 
earized Lorenz flow near EQ t equilibrium, see fig- 
ure 3.4(a). The unstable eigenplane of EQi is 
spanned by Ree (1) and Ime (1) . The stable eigen- 
vector e (3) . (b) Lorenz flow near the EQo equi- 
librium: unstable eigenvector e (1> , stable eigen- 



vectors e 
10" 



,(2) P (3) 



10 



-12 



. Trajectories initiated at distances 
10~ 13 away from the z-axis exit fi- 

d) p(2)l 



nite distance from EQo along the (e l ,e ( ') eigen- 
vectors plane. Due to the strong A {1) expansion, the 
EQ equilibrium is, for all practical purposes, un- 
reachable, and the EQ\ — > EQo heteroclinic con- 
nection never observed in simulations such as fig- 
ure 2.5. (E. Siminos; continued in figure 11.8.) 




(b) 



As all numerical plots of the Lorenz flow are here carried out for the Lorenz 
parameter choice cr = 10, b = 8/3, p = 28 , we note the values of these eigenvalues for 
future reference, 

EQ Q : (A m ,A {2 \A {i) ) = ( 11.83,- 2.666, -22.83) , . 

£Qi : (yu (1) + ;w (1) ,i (3) ) = ( 0.094 ±i 10.19, -13.85), ( ' 

as well as the rotation period T EQl = 2n/cj (1) about EQ\, and the associated expan- 
sion/contraction multipliers A® = exp(fi U) T EQl ) per a spiral-out turn: 

T EQl =0.6163, (A (1) ,A (3) ) = ( 1.060,1. 957 xl0~ 4 ). (4.40) 



We learn that the typical turnover time scale in this problem is of order T « T EQl « 1 
(and not, let us say, 1000, or 10~ 2 ). Combined with the contraction rate (4.41), this tells 
us that the Lorenz flow strongly contracts state space volumes, by factor of ~ 10~ 4 per 
mean turnover time. 

In the EQi neighborhood the unstable manifold trajectories slowly spiral out, 
with very small radial per-turn expansion multiplier A (1) 1.06, and very strong con- 
traction multiplier A (3) ^ 10~ 4 onto the unstable manifold, figure 4.6(a). This contraction 
confines, for all practical purposes, the Lorenz attractor to a 2-dimensional surface evi- 
dent in the section figure 3.4. 

In the x EQo = (0, 0, 0) equilibrium neighborhood the extremely strong A (3) = 
-23 contraction along the e (3) direction confines the hyperbolic dynamics near EQo to 
the plane spanned by the unstable eigenvector e (1) , with A m 12, and the slowest 
contraction rate eigenvector e (2) along the z-axis, with A (2} = -3. In this plane the strong 
expansion along e (1) overwhelms the slow A {2) -3 contraction down the z-axis, making 
it extremely unlikely for a random trajectory to approach EQ , figure 4.6(b). Thus 
linearization suffices to describe analytically the singular dip in the Poincare sections 
of figure 3.4, and the empirical scarcity of trajectories close to EQo. (continued in 
example 4.8) 

(E. Siminos and J. Halcrow) 



Example 4.7 Lorenz flow: Global portrait. (continued from example 4.6) As the 
EQi unstable manifold spirals out, the strip that starts out in the section above EQ\ in 
figure 3.4 cuts across the z-axis invariant subspace. This strip necessarily contains a 
heteroclinic orbit that hits the z-axis head on, and in infinite time (but exponentially fast) 
descends all the way to EQo. 
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How? As in the neighborhood of the EQ equilibrium the dynamics is linear 
(see figure 4.6(a)), there is no need to integrate numerically the final segment of the 
heteroclinic connection - it is sufficient to bring a trajectory a small distance away from 
EQo, continue analytically to a small distance beyond EQo, then resume the numerical 
integration. 

What happens next? Trajectories to the left of z-axis shoot off along the e (1) 
direction, and those to the right along -e (I) . As along the e (1) direction xy > 0, the 
nonlinear term in the z equation (2.13) bends both branches of the EQo unstable man- 
ifold W"(EQ ) upwards. Then ... - never mind. Best to postpone the completion of 
this narrative to example 9.14, where the discrete symmetry of Lorenz flow will help us 
streamline the analysis. As we shall show, what we already know about the 3 equilib- 
ria and their stable/unstable manifolds suffices to completely pin down the topology of 
Lorenz flow, (continued in example 9. 14) 

(E. Siminos and J. Halcrow) 



Example 4.8 Lorenz flow state space contraction: (continued from exam- 
ple 4.6) It follows from (4.30) and (4.29) that Lorenz flow is volume contracting, 

3 

divi = J- ( '\x, f) = -a-b-l, (4.41) 

!=1 

at a constant, coordinate- and p-independent rate, set by Lorenz to divi = -13.66 . As 
for periodic orbits and for long time averages there is no contraction/expansion along 
the flow, /l (ll) = 0, and the sum of A {l) is constant by (4.41), there is only one independent 
exponent A {>) to compute, (continued in example 4. 7) click to return: p. ?? 



Example 4.9 Stability of a 1 -dimensional map: Consider the orbit {. . .,jc-i,xo,jci,jc2, ■ • •) 

of a 1 -dimensional map x„ +l = f(x„). Since point x„ is carried into point x„ +l , in study- 
ing linear stability (and higher derivatives) of the map it is often convenient to deploy 
a local coordinate systems z„ centered on the orbit points x„, together with a notation 
for the map, its derivative, and, by the chain rule, the derivative of the kth iterate f k 
evaluated at the point x a , 

X = X a + Za , fa(Za) = f(x„ + Za) 
fa = f'tXa) 

A(x ,k) = f a ' =f a+k _ l ■■■fUJ'a, k>2. (4.42) 



stability - 17apr2013 



ChaosBook.org versionl4.4.1, Apr 21 2013 



EXERCISES 



96 



Here a is the label of point x a , and the label a+l is a shorthand for the next point b on 
the orbit of x a , xb = x„+\ = f(x a ). For example, a period-3 periodic point in figure 4.7 
might have label a = Oil, and by xno = /Oon) the next point labelis b - 110. click to return: p. ?? 



Example 4.10 Henon map Jacobian matrix: For the Henon map (3. 1 7) the Jaco- 
bian matrix for the nth iterate of the map is 

M n (x ) = f] ( ~ 2a 1 Xm b Q ), x m = /TOtcyo) - (4.43) 

The determinant of the Henon one time-step Jacobian matrix (4.43) is constant, 

detM = A[A 2 = (4.44) 

so in this case only one eigenvalue Ai = -b/A 2 needs to be determined. This is not 
an accident; a constant Jacobian was one of desiderata that led Henon to construct a 
map of this particular form. click to return: p. ?? 



Exercises 



4. 1 . Trace-log of a matrix. Prove that 



detM = e trlnM . 

for an arbitrary nonsingular finite dimensional matrix M, 
detM + 0. 

4.2. Stability, diagonal case. Verify the relation (4.32) 

/ = e tA = try Ao u , A D = UAIT 1 . 

4.3. State space volume contraction. 

(a) Compute the Rossler flow volume contraction rate 
at the equilibria. 

(b) Study numerically the instantaneous d,v, along a 
typical trajectory on the Rossler attractor; color- 
code the points on the trajectory by the sign (and 
perhaps the magnitude) of 3,v,-. If you see regions 
of local expansion, explain them. 

(c) (optional) color-code the points on the trajectory 
by the sign (and perhaps the magnitude) of <5,v; - 

(d) Compute numerically the average contraction rate 
(4.29) along a typical trajectory on the Rossler at- 
tractor. Plot it as a function of time. 



(e) Argue on basis of your results that this attractor is 
of dimension smaller than the state space d — 3. 

(f) (optional) Start some trajectories on the escape 
side of the outer equilibrium, color-code the points 
on the trajectory. Is the flow volume contracting? 

(continued in exercise 20.1 1) 

4.4. Topology of the Rossler flow. (continuation of exer- 
cise 3.1) 

(a) Show that equation |det(A - Al)\ = for Rossler 
flow in the notation of exercise 2.8 can be written 



A 3 +A 2 c( P T -e)+A(p ± /e+l-c 2 e P T )+cyfD = 

(b) Solve (4.45) for eigenvalues A ± for each equilib- 
rium as an expansion in powers of e. Derive 

/tj" = -c + ec/{c 2 + 1) + o(e) 



A- = ec 3 /[2(c 2 + l)] + o(e 2 ) 
6 2 = l + e/[2(c 2 + l)] + o(e) 
A\ =ce(l -e) + o(e 3 ) 
4 = - e 5 c 2 /2 + o(e 6 ) 
6»+ = VI + l/e(l +o(e)) 



(4.46) 
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Compare with exact eigenvalues. What are dy- 
namical implications of the extravagant value of 
/Lj"? (continued as exercise 13.10) 

(R. Paskauskas) 

4.5. Time-ordered exponentials. Given a time dependent 
matrix V(t) check that the time-ordered exponential 

<W(f) = TeX'^« 

may be written as 



with the initial condition 1/(0) = 1 . 

4.6. A contracting baker's map. Consider a contracting 
(or 'dissipative') baker's map, acting on a unit square 
[0, 1] 2 = [0,1] x [0,1], defined by 



rrS Jo jo 



f dh--- f 

Jo Jo 



dt m V(h)---V{t m ) 



and verify, by using this representation, that satis- 
fies the equation 

<ii(t) = V(t)ti(t), 



y n +\ 

x n +\ 
y n +\ 



x n /3 
2y„ 



y„ < 1/2 



*„/3 + l/2 
2y„ - 1 



y n > 1/2. 



This map shrinks strips by a factor of 1/3 in the x- 
direction, and then stretches (and folds) them by a factor 
of 2 in the y-direction. 

By how much does the state space volume contract for 
one iteration of the map? 
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